// http://www.geometer.org/puzzles/soma.c

// Copyright 2005 Thomas R. Davis
//
// This file is part of Soma.
//
// Soma is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Soma is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Soma; if not, write to the Free Software Foundation,
// Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

#include <stdio.h>
#include <stdlib.h>

// This code solves the soma cube by fitting the 7 pieces into
// a 3x3x3 block.  It is done purely by brute force.

// Note that this, in a sense, double-counts
// the solutions, since every piece is a mirror
// image of itself, except that C and D are
// mirror images of each other.  If the pieces
// are painted different colors, there are 480
// solutions; if the same, there are only 240,
// since you can mirror any solution and the D
// and C cubes swap jobs.

// Following are the coordinates for one of
// each of the typical blocks.  The only thing
// that's important here is that the block
// have one cube with coordinates (0,0,0).

int Apiece[4][3] = {
{0, 0, 0}, {0, 0, 0}, {0, 1, 0}, {1, 0, 0}
};

int Bpiece[4][3] = {
{0, 0, 0}, {0, 0, 1}, {0, 1, 0}, {1, 0, 0}
};

int Cpiece[4][3] = {
{0, 0, 0}, {0, 1, 1}, {0, 1, 0}, {1, 0, 0}
};

int Dpiece[4][3] = {
{0, 0, 0}, {1, 0, 1}, {0, 1, 0}, {1, 0, 0}
};

int Epiece[4][3] = {
{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {1, 1, 0}
};

int Fpiece[4][3] = {
{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}
};

int Gpiece[4][3] = {
{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 2, 0}
};

// These arrays hold the encoded positions of
// the cubes in a bitmapped way.  Bits 0 through
// 26 correspond to the locations in a 3x3x3
// cube.  They are large since all configurations
// are generated, but then sorted and duplicates
// are erased later.

int Aposns[500];
int Bposns[500];
int Cposns[500];
int Dposns[500];
int Eposns[500];
int Fposns[500];
int Gposns[500];

// counts of piece configurations:

int Ac, Bc, Cc, Dc, Ec, Fc, Gc;

// These are all the matrices for 3x3x3
// rotations in space generated by the
// 90 degree turns.  But not reflections.

#define M -1

int m[24][3][3] = {
{{1,0,0},{0,1,0},{0,0,1}},
{{0,0,1},{0,1,0},{M,0,0}},
{{0,0,1},{0,M,0},{1,0,0}},
{{0,0,1},{1,0,0},{0,1,0}},
{{0,0,1},{M,0,0},{0,M,0}},
{{0,0,M},{0,1,0},{1,0,0}},
{{0,0,M},{0,M,0},{M,0,0}},
{{0,0,M},{1,0,0},{0,M,0}},
{{0,0,M},{M,0,0},{0,1,0}},
{{0,1,0},{0,0,1},{1,0,0}},
{{0,1,0},{0,0,M},{M,0,0}},
{{0,1,0},{1,0,0},{0,0,M}},
{{0,1,0},{M,0,0},{0,0,1}},
{{0,M,0},{0,0,1},{M,0,0}},
{{0,M,0},{0,0,M},{1,0,0}},
{{0,M,0},{1,0,0},{0,0,1}},
{{0,M,0},{M,0,0},{0,0,M}},
{{1,0,0},{0,0,1},{0,M,0}},
{{1,0,0},{0,0,M},{0,1,0}},
{{1,0,0},{0,M,0},{0,0,M}},
{{M,0,0},{0,0,1},{0,1,0}},
{{M,0,0},{0,0,M},{0,M,0}},
{{M,0,0},{0,1,0},{0,0,M}},
{{M,0,0},{0,M,0},{0,0,1}},
};

// for debugging:
void sanitycheckrots()
{
  int i, det;
  
  for (i = 0; i < 24; i++) {
    det =  m[i][0][0]*m[i][1][1]*m[i][2][2]
         + m[i][1][0]*m[i][2][1]*m[i][0][2]
         + m[i][2][0]*m[i][0][1]*m[i][1][2]
         - m[i][0][0]*m[i][2][1]*m[i][1][2]
         - m[i][1][0]*m[i][0][1]*m[i][2][2]
         - m[i][2][0]*m[i][1][1]*m[i][0][2];
    if (det != 1)
        fprintf(stderr, "%d, error!\n", i);
  }
}

void rotate(int in[4][3], int out[4][3], int rot)
{
  int i, j, k;
  
  for (i = 0; i < 4; i++) {
    for (j = 0; j < 3; j++) {
      out[i][j] = 0;
      for (k = 0; k < 3; k++)
        out[i][j] += in[i][k]*m[rot][k][j];
    }
  }
}

// mix: finds a unique number from 0-26 for
// a 3-dimensional set of coordinates.

int mix(int a, int b, int c)
{
  return a*9 + b*3 + c;
}

// compares two numbers for qsort

int comp(const void *a, const void *b)
{
  if (*(int *)a < *(int *)b) return -1;
  if (*(int *)a == *(int *)b) return 0;
  return 1;
}

// generates all rotations and translations
// of a piece.  Then it sorts and eliminates
// the duplicates.

int fillslots(int piece[4][3], int outpiece[500])
{
  int count = 0;
  int i, j, k, l, r, ok, last, ind;
  int out[4][3];
  
  for (r = 0; r < 24; r++)
    for (i = 0; i < 3; i++)
      for (j = 0; j < 3; j++)
        for (k = 0; k < 3; k++) {
          rotate(piece, out, r);
          ok = 1;
          for (l = 0; l < 4; l++) {
            out[l][0] += i;
            out[l][1] += j;
            out[l][2] += k;
            if (out[l][0]<0 || out[l][0]>2) ok=0;
            if (out[l][1]<0 || out[l][1]>2) ok=0;
            if (out[l][2]<0 || out[l][2]>2) ok=0;
          }
          if (ok) {
            outpiece[count] = 0;
            for (l = 0; l < 4; l++)
              outpiece[count] |= (1 << 
                mix(out[l][0], out[l][1], out[l][2]));
            count++;
          }
        }
  qsort(outpiece, count, sizeof(int), comp);
  
  // nuke the dups:
  
  last = -1;
  ind = 0;

  for (i = 0; i < count; i++)
    if (outpiece[i] != last) {
      outpiece[ind++] = last = outpiece[i];
    }
  return ind;
}

int pos[4][3] = {
{0, 0, 0},
{1, 0, 0},
{1, 1, 0},
{1, 1, 1}
};

// Puts one particular piece, the C piece, in
// the only 4 possible positions that are unique
// up to rotation.  There's no need to treat C
// differently, except that this eliminates all
// of the symmetrical solutions.  Or all but half
// of them, depending on whether you think the
// mirror images that swap C and D are the same.

int fillspecial(int piece[4][3], int outpiece[500])
{
  int i, j, k, l, count = 0;
  int out[4][3];
  
  for (i = 0; i < 4; i++) {
          rotate(piece, out, 0);
          for (l = 0; l < 4; l++) {
            out[l][0] += pos[i][0];
            out[l][1] += pos[i][1];
            out[l][2] += pos[i][2];
          }
          outpiece[count] = 0;
          for (l = 0; l < 4; l++)
            outpiece[count] |= (1 << 
                mix(out[l][0], out[l][1], out[l][2]));
          count++;
      }
   return count;
}

int Cube = 0;

// checks to see if a piece fits; if so, it puts it
// in and returns 1; otherwise, it returns 0.

int load(int a)
{
  if (a & Cube) return 0;
  Cube = Cube | a;
  return 1;
}

// unload nukes a piece from the cube

void unload(int a)
{
  Cube = Cube & (~a);
}

// dumppos dumps a found configuration

dumppos(int a, int b, int c, int d, int e, int f, int g)
{
  int i, probe;
  
  for (i = 0; i < 27; i++) {
    probe = 1 << i;
    if (probe & Aposns[a]) printf("a");
    if (probe & Bposns[b]) printf("b");
    if (probe & Cposns[c]) printf("c");
    if (probe & Dposns[d]) printf("d");
    if (probe & Eposns[e]) printf("e");
    if (probe & Fposns[f]) printf("f");
    if (probe & Gposns[g]) printf("g");
    if (i%3 == 2) printf(" ");
    if (i%9 == 8) printf("\n");
  }
  printf("\n");
}

// main looks for all the solutions.  The order
// of the nesting was chosen to have the largest
// sets of possible positions in the inner loops.
// For the soma cube this isn't too important,
// but it would be for larger puzzles and a
// brute-force solution such as this.

int main()
{
  // sanitycheckrots();
  int a, b, c, d, e, f, g;
  int solutions = 0;
  
  Ac = fillslots(Apiece, Aposns);
  Bc = fillslots(Bpiece, Bposns);
  Cc = fillspecial(Cpiece, Cposns);
  Dc = fillslots(Dpiece, Dposns);
  Ec = fillslots(Epiece, Eposns);
  Fc = fillslots(Fpiece, Fposns);
  Gc = fillslots(Gpiece, Gposns);
  
  //printf("%d %d %d %d %d %d %d\n",
    //Ac, Bc, Cc, Dc, Ec, Fc, Gc);
    
  //printf("%x %x %x %x %x %x %x\n",
    //Aposns[0], Bposns[0], Cposns[0],
    //Dposns[0], Eposns[0], Fposns[0], Gposns[0]);
    
  for (c = 0; c < Cc; c++) {
    load(Cposns[c]);
    for (b = 0; b < Bc; b++) {
      if (0 == load(Bposns[b])) goto labb;
      for (f = 0; f < Fc; f++) {
        if (0 == load(Fposns[f])) goto labf;
        for (d = 0; d < Dc; d++) {
          if (0 == load(Dposns[d])) goto labd;
          for (e = 0; e < Ec; e++) {
            if (0 == load(Eposns[e])) goto labe;
            for (g = 0; g < Gc; g++) {
              if (0 == load(Gposns[g])) goto labg;
              for (a = 0; a < Ac; a++) {
                if (0 == load(Aposns[a])) goto laba;
                dumppos(a, b, c, d, e, f, g);
                solutions++;
                unload(Aposns[a]);
              laba:;}
              unload(Gposns[g]);
            labg:;}
            unload(Eposns[e]);
          labe:;}
          unload(Dposns[d]);
        labd:;}
        unload(Fposns[f]);
      labf:;}
      unload(Bposns[b]);
    labb:;}
    unload(Cposns[c]);
  }
  printf("Total solutions found: %d\n", solutions);
}
         
